#include "share/Usrdef.h"
#include "core/globalvar.h"
#include "core/Noise.h"
#include "share/parameters.h"
#include "share/gsl_app.h"
#include "share/MatAlgebra.h"

extern double		NOISE_AMPLITUDE;
extern double		CORRELATION_TIME;

double		* POWER_LIST;

///////////////////////////////////////////////////////////////////////////////////

void noise_ensemble::initialize(double minT, double maxT, int nt, double minW, double maxW, int nW, int num_path, int num_interpolation)
{
	cout << "begin initialization..." << endl;

	int i, j, idx;

	mintime = minT;
	maxtime = maxT;
	ntime 	= nt;
	dt		= (maxT - minT) / nt;

	minfrequency = minW;
	maxfrequency = maxW;
	nfrequency = nW;
	df		= (maxW - minW) / nW;

	npath 	= num_path;
	ninterpolation = num_interpolation;

	harmonic_mat_sin = new double [nfrequency * ntime];
	harmonic_mat_cos = new double [nfrequency * ntime];
	harmonic_mat_sin_interpolation = new double [nfrequency * ninterpolation];
	harmonic_mat_cos_interpolation = new double [nfrequency * ninterpolation];

	interpolation_time_list = new double [ninterpolation];
	bt_interpolation = new double [ninterpolation];

	double * wt_list = new double [nfrequency * ninterpolation];
	double inter_dt;

	double ti, wj;
	for( i = 0; i < ntime; i ++)
	{
		ti = i * dt;
		for(j = 0; j < nfrequency; j ++)
		{
			wj = j * df;
			idx = i * nfrequency + j;
			harmonic_mat_sin[ idx ] = sin(wj * ti);
			harmonic_mat_cos[ idx ] = cos(wj * ti);

		}
	}

	inter_dt = ( maxtime - mintime ) / ninterpolation;
	for( i = 0; i < ninterpolation; i ++)
	{
		ti = i * inter_dt;
		interpolation_time_list[i] = ti;
		for(j = 0; j < nfrequency; j ++)
		{
			wj = j * df;
			idx = i * nfrequency + j;
			wt_list[idx] = wj * ti;
		}
	}
	vdCos(nfrequency * ninterpolation, wt_list, harmonic_mat_cos_interpolation);
	vdSin(nfrequency * ninterpolation, wt_list, harmonic_mat_sin_interpolation);

	cout << "path initialized." << endl;

	delete [] wt_list;
}

void noise_ensemble::generate_bath(int s)
{
		const gsl_rng_type * T;
		gsl_rng * r;

		gsl_rng_env_setup();

		T = gsl_rng_default;
		r = gsl_rng_alloc (T);

		gsl_rng_set(r, s);
//////////////////////////////////////////////////////////////////////////
		int i, j, idx;

		amplitude_mat_alpha = new double [npath * nfrequency];
		amplitude_mat_beta = new double  [npath * nfrequency];

//		for(i = 0; i < npath; i ++)
//		{
//			for(j = 0; j < nfrequency; j ++)
//			{
//				idx = j * npath + i;
//				amplitude_mat_alpha[idx] = gsl_ran_gaussian( r, POWER_LIST[j] );
//				amplitude_mat_beta[idx] = gsl_ran_gaussian( r, POWER_LIST[j] );
//			}
//		}
		for(i = 0; i < npath; i ++)
		{
			for(j = 0; j < nfrequency; j ++)
			{
				idx = i * nfrequency + j;
				amplitude_mat_alpha[idx] = gsl_ran_gaussian( r, POWER_LIST[j] );
				amplitude_mat_beta[idx] = gsl_ran_gaussian( r, POWER_LIST[j] );
			}
		}

		gsl_rng_free (r);

//	    char amp_file[200];
//
//	    strcpy(amp_file, WORKSHOP);
//	    strcat(amp_file, DEBUG_PATH);
//	    strcat(amp_file, "AMP_");
//	    strcat(amp_file, PARA_CHAR);
//	    strcat(amp_file, ".dat");
//
//	    ofstream amp(amp_file);
//	    if(!amp) assert(0);
//	    idx = 0;
//		for(i = 0; i < npath; i ++)
//		{
//			for(j = 0; j < nfrequency; j ++)
//			{
//				amp << setw(12) << setiosflags(ios::scientific) << setprecision(10) << amplitude_mat_alpha[idx] <<"\t" << amplitude_mat_beta[idx] << endl;
//				idx ++;
//			}
//		}
//		amp.close();
}

void noise_ensemble::generate_interpolation_point(int path_idx, int power)
{
	int i, j;

	int m, n, lda, inc;
	double alpha, beta;
	m = nfrequency;
	n = ninterpolation;
	lda = m;
	inc = 1;
	alpha = CONST_ONE;
	beta = CONST_ZERO;
	dgemv(TRANS, &m, &n, &alpha, harmonic_mat_cos_interpolation, &lda, amplitude_mat_alpha + path_idx * nfrequency,
			&inc, &beta, bt_interpolation, &inc);

	beta = CONST_ONE;
	dgemv(TRANS, &m, &n, &alpha, harmonic_mat_sin_interpolation, &lda, amplitude_mat_beta + path_idx * nfrequency,
			&inc, &beta, bt_interpolation, &inc);

	double Z;
	if(power == 2)
	{
		for(i = 0; i < ninterpolation; i ++)
			bt_interpolation[i] = (bt_interpolation[i] * bt_interpolation[i]);
	}
	else if(power == 3)
	{
		for(i = 0; i < ninterpolation; i ++)
		{
			Z = OMEGA_Z + bt_interpolation[i];
			bt_interpolation[i] = sqrt(Z * Z + OMEGA_V * OMEGA_V);
		}
	}
}

double noise_ensemble::eval(int n, double tval)
{
	double res, wk;
	int k, idx;

	res = 0.0;
	for(k = 0; k <nfrequency; k++)
	{
		wk = k * df;
		idx =  n * npath + k;
		res += amplitude_mat_alpha[idx] * cos(wk * tval) + amplitude_mat_beta[idx] * sin(wk * tval);
	}
	return res;
}

void noise_ensemble::calc_correlation_func(double * res)
{
	int i, inc = 1;
	bt_mat = new double [npath * ntime];

	calc_quantity(harmonic_mat_sin, harmonic_mat_cos, bt_mat);

	for(i = 0; i < ntime; i ++)
	{
		res[i] = ddot(&npath, bt_mat, &inc, bt_mat+i*npath, &inc) / npath;
	}

}

void noise_ensemble::set_phase_mat(double * ph)
{
	int i;
	phase_mat = new double [npath * ntime];
	for(i = 0; i < npath * ntime; i++)
		phase_mat[i] = ph[i];
}

void noise_ensemble::calc_dephasing(double * coherence_re, double * coherence_im)
{
	int i, j, idx;
	double ti;
	double * cos_phase_mat, * sin_phase_mat;

	cos_phase_mat = new double[NTIME * NPATH];
	sin_phase_mat = new double[NTIME * NPATH];

	vdCos(NTIME * NPATH, phase_mat, cos_phase_mat);
	vdSin(NTIME * NPATH, phase_mat, sin_phase_mat);
	for(i = 0; i < NTIME; i++)
	{
		ti = i * TIME_STEP;

		coherence_re[i] = 0.0;
		coherence_im[i] = 0.0;
		for(j = 0; j < NPATH; j++)
		{
			idx = j * NTIME + i;//i* NPATH + j;
			coherence_re[i] += cos_phase_mat[idx];
			coherence_im[i] += sin_phase_mat[idx];
		}
		coherence_re[i] /= NPATH;
		coherence_im[i] /= NPATH;
	}
	delete [] cos_phase_mat;
	delete [] sin_phase_mat;
}

void noise_ensemble::calc_quantity(double * sin_mat, double * cos_mat, double * res_mat)
{
	MKL_real_mat_mult(npath, ntime, nfrequency, TRANS, NO_TRANS, CONST_ONE, amplitude_mat_alpha, cos_mat, CONST_ZERO, res_mat);
	MKL_real_mat_mult(npath, ntime, nfrequency, TRANS, NO_TRANS, CONST_ONE, amplitude_mat_beta, sin_mat, CONST_ONE, res_mat);
}
